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In  use  at  the  Naval  Undersea  Warfare  Center,  San  Diego  Division. This 
note  is  not  to  he  considered  an  official  NUWC  report.  Its  purpose  is 
to  docvunent  the  prediction  technique  for  the  benefit  of  those  who  may 
have  need  to  use  it. 


The  work  described  in  this  technical  note  has  been  supported  under 
NAVSHIPS  Exploratory  Development  subproject  SF  101  03  21,  Task  8704, 
and  by  Independent  Research  funds  under  NAVSHIPS  subproject  ZR  011  01  01. 


A CONTINUOUS-SLOPE  CURVE  FITTING  ROUTINE 
FOR  USE  IN 

SONAR  PERFORt4ANCE  PREDICTION 

Ray  tracing  techniques  have  long  been  used  to  estimate  the  propaga- 
tion loss  of  sound  in  the  ocean.  The  availability  of  high-speed  digital 
computers  has  made  practical  the  application  of  such  techniques  to  in- 
creasingly complex  environmental  conditions.  The  "velocity  profile", 
which  specifies  the  dependence  of  sound  velocity  C on  depth  below  the 
surface  of  the  ocean  Z,  is  the  most  critical  environmental  input  to 
any  ray -tracing  routine.  The  velocity  profile  is  specified  as  a set  of 
points  input  to  a digital  computer  program.  A profile  for 

which  ray-path  solutions  can  be  obtained  in  closed  form  is  then  fit  to 
the  set  of  points.  A relatively  simple  solution  can  be  obtained  for  a 
profile  consisting  of  consecutively  joined  straight-line  segments  passing 

through  all  the  points  (C.,Z.).  However,  the  slope  discontinuities  which 
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inevitably  occur  in  such  striaght-line-segment  curves  lead  to  the  predic- 
tion of  non-physical  singularities  in  the  sound  intensity.  More  satisfac- 

t 

tory  predictions  have  been  obtained  by  Pedersen,  Gordon,  and  Keith, 
using  a continuous  slope  curve  to  fit  a given  velocity  profile.  The 
fitted  curve  consists  of  joined  segments,  each  of  which  satisfies  an 
equation  of  the  form 

a + bZ  + cZ^  = C-2. 

Pedersen,  Gordon,  and  Keith^  have  worked  out  the  ray-trace  theory  in 
detail  for  this  case  and  have  written  flow  charts  for  a computer  program 


to  map  sound  intensities.  The  entire  ray  trace  package  is  being  programmed 


In  ETRAN*  by  NUWC  Code  D556.  This  report  will  describe  the  curve-fitting 
part  of  the  overall  ray  trace  package.  The  remainder  of  the  ray  trace 
package  will  be  described  in  a subsequent  report. 

The  "computer  notation"  used  in  the  balance  of  this  report  will  treat 
as  independent  variables  quantities  that,  in  mathematical  notation,  would 
normally  be  written  as  subscripts,  e.g.,  and  will  be  written 

C(I)  and  Z(I). 

The  curve-fitting  routine  divides  the  ocean  between  the  shallowest 
point  and  the  deepest  point  of  the  input  velocity  profile  into  layers 
bounded  by  horizontal  planes,  as  indicated  in  Figure  1.  Depths  of  the 
bounding  planes  will  be  denoted  by  ZU.  The  layer  numbered  IL  extends 
from  ZU(IL)  downward  to  ZU(lL+l).  Parameters  AA,  BB,  and  CC  are 
computed  for  each  layer  such  that  the  condition 

|AA(IL)  + BB(IL)  X ZS  + CC(IL)  x ZS2-CS-2|<EM 

is  satisfied  for  each  point  (CS,ZS)  of  the  input  profile  for  which 
ZU(lL)sZSiZU(lL+l).  Here,  EM  is  an  input  parameter  specifying  the 
maximum  distance  in  the  (C“^,Z)  plane  that  the  fitted  curve  is  allowed 
to  deviate  from  any  point  of  the  input  profile.  The  coefficients  AA(IL), 
BB{IL),  CC(IL),  and  AA(lL+l),  BB(lL+l),  CC(lL+l)  must  satisfy  the 
condition  that  the  slope  of  the  segment  numbered  IL  equals  that  of  the 
segment  numbered  IL+1  at  the  junction  depth  ZU(lL+l). 

*ETRAN  is  a scientific  programming  language  developed  by  NUVv'C  Code  D556 
(formerly  NEL  Code  3110D)  used  in  Sonar  Performance  Prediction. 
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The  primary  task  of  the  curve-fitting  program,  then,  is  to  determine 
layer  interface  depths  ZU  and  the  coefficients  AA,  BB,  and  CC  for  each 
layer.  The  first  step  is  to  make  a least-squares  fit  of  the  equation 

a + hZ  + cZ^  = C"^ 

to  the  first  four  points  of  the  velocity  profile. 

If  the  condition 

la  + b X ZS  + c X ZS2  - CS-^ 1 <EM  (l) 

is  not  satisfied  for  each  of  the  four  points,  then  an  exact  fit  is  made  to 
the  first  three  points.  If  the  four-point  fit  is  satisfactory,  then  a least- 
squares  fit  is  made  to  the  first  five  points  and  each  point  is  tested  against 
equation  (l).  The  nmber  of  points  is  increased  by  one  until  an  unsatisfactory 
fit  is  obtained,  or  until  all  points  have  been  fitted  satisfactorily.  The 
coefficients  of  the  last  last  satisfactory  fit  are  stored  as  AA(0),  BB(0), 

CC(0).  The  depth  of  the  deepest  point  used  in  the  last  satisfactory  fit  is 
stored  as  ZU(3).  The  depth  of  the  next  shallower  point  is  stored  as  ZU(l). 
Layer  0 thus  extends  from  the  depth  of  the  shallowest  point,  ZU(0),  to  a 
depth  ZU(l).  Layers  1 and  2 are  contained  between  depths  ZU{l)  and  ZU(3). 

The  program  returns  later  to  compute  ZU(2)  and  the  coefficients  AA,  BB, 
and  CC  for  layers  1 and  2. 

Starting  with  the  point  at  depth  ZU(3),  the  program  fits  another  group 
of  points  using  the  same  procedure  that  was  used  to  fit  the  first  group.  .The 
quantities  ZU(1*),  ZU(6),  AA(3),  BB(3),  and  CC(3)  are  computed. 
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This  procedure  is  continued  until  the  deepest  point  has  been  fit.  A 
special  situation  arises  if  a fit  including  the  next-to-deepest  point  is 
satisfactory,  but  a fit  including  the  deepest  point  as  well  is  unsatisfactory. 

In  this  case,  a curve  is  constructed  which  passes  through  the  deepest  of  the 
input  points  and  which  meets  the  preceding  curve  with  matching  slope  at  a depth 
equal  that  of  the  next-to-deepest  point. 

After  all  the  points  have  been  fitted,  the  program  goes  back  and  fills  in 
all  the  gaps  with  bridging  curves.  The  first  bridging  curve  extends  through 
layers  1 and  2.  Recall  that  ZU(l),  ZU(3),  and  the  coefficients  AA,  BB, 
and  CC  for  layers  0 and  3 have  already  been  determined.  It  is  necessary  to 
compute  ZU(2)  and  AA,  BB,  and  CC  for  layers  1 and  2.  The  depth  ZU(2) 
is  taken  as  the  depth  of  intersection  of  the  tangent  to  curve  0 at  depth  ZU(l) 
and  the  tangent  to  curve  3 at  depth  ZU(3),  provided  that  this  depth  of  inter- 
section occurs  at  a depth  greater  than  ZU(3),  then  ZLJ(2)  = 1/3  ZU(l)  + 2/3  ZU(3). 
Coefficients  are  now  computed  for  the  bridging  curves  1 and  2.  The  conditions 
which  must  be  satisfied  are  that  curves  0 and  1 meet  and  have  equal  slopes  at 
depth  ZU(l),  that  curves  2 and  3 meet  and  have  equal  slopes  at  depth  ZU(3), 
and  that  curves  1 and  2 meet  with  equal  slopes  at  depth  ZU(2). 

Other  necessary  bridging  curves  are  constructed  in  the  same  manner. 

Using  values  determined  for  ZU,  AA,  BB,  and  CC,  additional  parameters 
P2,  P3,  p4,  and  P5,  which  are  required  as  inputs  to  the  ray  tracing  program, 

are  computed  for  each  layer.  If  a positive-gradient  surface  channel  is  present, 
the  channel  depth  and  average  velocity  gradient  in  the  channel  are  determined. 
Finally,  points  on  the  fitted  velocity  profile  are  printed  out  at  closely- 
spaced  depth  intervals. 
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A complete  list  of  inputs  and  the  units  of  each  is: 

CS  (velocity)  , 

ZS  (length)  / ' 

NM 

EM  (velocity-^) 

CH  (length) , 

where 

CS  and  ZS  are  arrays  giving  the  points  of  the  input  velocity  of 
profile  in  order  of  increasing  depth. 

NM  is  the  number  of  points  in  the  arrays  ZS  and  CS. 

EM  is  the  maximum  distance  in  the  (C~^,Z)  plane  that  the  fitted  curve 
is  allowed  to  deviate  from  any  of  the  given  points.  It  is  desirable 
to  determine  EM  in  terms  of  deviation  allowed  in  the  (C,Z)  plane. 

Note  that 

d(c-2)  = _ . 

c^ 

AC 

Let  C^  be  a point  in  the  velocity  profile.  If  — <<  1,  then 

o 

ls(c-2)|  ; I I 

C 2 

O 

Letting  AC  - aximum  deviation  allowed  in  the  (C,Z)  plane, 


EM 


2 X AC 


o 


If  the  units  for  C are  ft/sec,  and  we  take  C^  = 1<950  ft/sec,  then 

EM  = 1.6t  X 10-*^ AC. 
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If  the  units  are  yds/sec 


EM  = U.5  X 10“^^ AC.  , • 

Pedersen,  Gordon,  and  Keith^  have  suggested  the  value  .2  yds/sec,  or  .6  ft/sec, 
for  AC.  These  numbers  give 

EM  = 1.00  X 10“^^  (units  of  C = ft/sec) 

EM  = 9.0  X 10-^^  (units  of  C = yds/sec). 

CH  is  the  depth  interval  at  which  points  from  the  fitted  profile  are 
printed  out.  Typically  specified  values  are  10  ft/sec  or  5 yds/sec. 

The  outputs  are: 


ILM 

ZP 

IL 

CP 

ZUI 

CL 

Cl 

ZL 

P2I 

GRD 

P3I 

VA 

pUi 

VB 

P5I 

VC 

"Profile  items"  (one  value  for  each  profile)  are: 

ILM  - number  of  layers  in  profile 
CL  - absolute  maximum  velocity  of  fitted  profile 
ZL  - depth  at  which  velocity  is  absolute  maximum 
GRD  - average  velocity  gradient  in  surface  channel. 


6 


"Layer  items"  (one  value  for  each  layer)  are: 

IL  = layer  number 

VA  = coefficient  AA  ^ . 

VB  = coefficient  BB 
VC  = coefficient  CC 

ZUI  = depth  of  upper  boundary  of  layer 
Cl  = velocity  at  upper  boundary  of.  layer 
P2I 
P3I 

' special  parameters  for  input  to  ray  tracing  program 

pUi 

P51 

"Profile  print-out  items" 

ZP 

points  taken  from  fitted  profile  at  depth  intervals  of  CH 
CP 

The  dotted  line  in  Figure  2A  represents  a curve  fitted  to  a oi-t  of  points 
denoted  by  X's.  The  points  are: 


Point  No. 

ZS(ft) 

CS(ft/sec) 

0 

0 

l)9‘+5.5 

1 

32.8 

1)91)6.1 

2 

65 . 6 

1)91)6.6 

3 

98.k 

l)9i)7.1 

h 

l6h 

1)951.9 

5 

21)6 

1)951.5 

6 

328 

1)989.8 

7 

1)92 

1)905.8 

8 

656 

1)888.0 
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Figure  2 


A least  squares  fit  to  the  set  of  points  (0,1, 2, 3)  was  first  made.  This 
fit  was  satisfactory  (with  EM  = 1.00  x 10“^^),  so  a fit  was  made  to  the  set 
(0,1 ,2 ,3 ,lt ) . This  fit  was  also  satisfactory.  A fit  . to  (0,1, 2, 3, ^,5)  was 
xinsatisfactory . The  values  ZU(o)  = 0,  ZU(l)  = 98.^,  and  ZU(3)  = 16!+  were 
stored.  Coefficients  determined  for  the  fit  to  (0,1,2,3,^<)  were  stored  as 
AA(0),  BB(0),  CC(o).  a four-point  fit  was  made  to  (U,5,6,T).  This  fit  was 
unsatisfactory,  so  an  exact  fit  was  made  to  (^,5,6).  The  values  ZU(U)  = 2U6 
and  ZU(6)  = 328  were  stored.  A three-point  fit  was  made  to  (6,7,8),  the 
last  three  input  points.  Next,  the  Bridging  curves  were  computed.  The 
quantity  ZU(2)  = 11*2.1  was  determined;  then,  the  coeff-cients  AA,  BB,  CC 
for  layers  1 and  2 were  computed.  Similarly  ZU(5)  = 300.7  was  determined; 
then,  AA,  BB,  and  CC  for  layers  1*  and  5 were  computed. 

The  X's  in  Figure  2A  represent  points  taken  at  standard  Nansen  depths 
from  a continuous  velocity  profile  determined  from  oceanographic  data.  Figure 
2B  compares  this  profile,  represented  by  a solid  line,  with  the  fitted  curve 
shown  in  Figure  2A,  represented  by  a dotted  line.  The  fact  that  there  are 
noticeable  differences  of  up  to  3 ft/sec  in  the  two  curves  is  due  to  the  large 
depth  intervals  between  the  points  to  which  the  dotted  curve  v/as  fitted.  This 
illustrates  the  importance  of  using  bathythermograph  or  other  continuously 
sampled  data  to  supplement  Nansen  Cast  data,  since  Nansen  sampling  intervals 
are  normally  too  gross  to  be  used  in  reproducing  profile  fine  structure.  A 
curve  fit  to  a set  of  points  taken  from  the  original  profile  at  depth  intervals 
of  25'  did  not  deviate  from  the  original  profile  by  more  than  0.6  ft /sec.  It 
is  evident  that,  if,  the  depth  interval  between  points  is  taken  to  be  small 
enough,  the  fitted  curve  can  be  made  to  correspond  closely  the  true  velocity 
profile  of  the  ocean. 
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